library(rstan)

n <- c(129, 35, 228, 84, 291, 270, 46, 293, 241, 105, 353, 250, 41)
y <- c(4,1,18,7,24,16,6,19,15,13,25,11,4)
I <- length(n)
eta_star <- 1.0

data <- list(
  n = n,
  y = y,
  I = I,
  eta_star = eta_star
)


fit <- stan(  file = 'hs_r_stan/hs_r.stan',
  data = data,
  iter = 10000,
  chains = 4,
  seed = 667,
  warmup = 1000
)

print(fit, pars = c("mu", "eta", "p", "y_pred"))


posterior <- as.data.frame(fit)

par(mfrow = c(2, 2))
hist(posterior$mu, main="Posterior of mu", xlab="mu", breaks=30, col="lightblue")
hist(posterior$eta, main="Posterior of eta", xlab="eta", breaks=30, col="lightgreen")
for(i in 1:I){
  hist(posterior[[paste0("p[", i, "]")]], 
       main=paste("Posterior of p[", i, "]"), 
       xlab=paste("p[", i, "]"), 
       breaks=30, col="lightyellow",
       y_lim=c(0,50))
}

y_pred_mean <- colMeans(posterior[, grep("^y_pred", names(posterior))])
plot(y, y_pred_mean, 
     main="Observed vs Predicted", 
     xlab="Observed y", 
     ylab="Predicted y", 
     pch=19, col="blue")

abline(0, 1, col="red", lwd=2)

par(mfrow=c(1,1))


253103.74+188754.78


